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We apply non-local density functional formalism to describe an equilibrium distribution of the 
water-like fluid in the asymmetric nanoscale junction presenting an atomic force microscope (AFM) 
tip dwelling above an arbitrary surface. The hydrogen bonding dominating in intermolecular at- 
traction is modeled as a square well potential with two adjustable parameters (energy and length) 
characterizing well's depth and width. A water meniscus formed inside nanoscale junction is explic- 
itly described for different humidity. Furthermore, we suggest a simple approach using polymolecular 
adsorption isotherms for the evaluation of an energetic parameter characterizing fluid (water) attrac- 
tion to substrate. This model can be easily generalized for more complex geometries and effective 
intermolecular potentials. Our study establishes a framework for the density functional description 
of fluid with orientational anisotropy induced by the non-uniform external electric field. 

PACS numbers: 



I. INTRODUCTION 

An equilibrium behavior of a fluid near geometrically 
nonuniform solid surfaces exhibits various peculiarities 
related to confinement and spatially varying external po- 
tentials. A quantitative understanding of water conden- 
sation phenomena in proximity of nano-asperities un- 
der ambient humidity is important for different areas 
of research including scanning probe microscopy (SPM), 
nano-patterning, adhesion and friction at macro- and 
nanoscale. Specifically, the formation of liquid films and 
bridges near an asperity of sub-micro- or nano-scale cur- 
vature needs to be adequately described. 

An importance of the water condensation phenomenon 
from ambient atmosphere into nanoscale junction formed 
by AFM tip dwelling above surface has been recognized 
in a number of works The formation of 
the water meniscus affects strongly the force-distance 
measurements and AFM imaging resolutions^ provides 
nanoscale electrochemical cell for SPM oxidation of semi- 
conductors and metals)^ serves as ink transport chan- 
nel for dip-pen nanolithography^ and plays an essential 
role in AFM-assisted electrostatic nanolithography (AF- 
MEN) of thin polymeric films/iVi 10 

Direct experimental observations of meniscus forma- 
tion in the AFM tip - surface junction are complicated 
due to nanometer-scale size of the region of interest. 
There are no optical tools available for this purpose due 
to the wavelength limitation. Indirect estimations of 
the meniscus's size can be made using SPM oxidation, 3 
noncontact AFM imaging 1 ^^ and substrate's dissolu- 
tion in water^ Several theoretical estimations of the 
meniscus size and shape based on the macroscopic Kelvin 
equationi^i^ and molecular level grand canonical Monte 
Carlo simulationsiSiiLiS have been performed. A macro- 
scopic phenomenological approach based on the modifi- 
cations of the Kelvin equation has been suggested to de- 
scribe water attracted to biased AFM tip>i£ The model 
misses molecular level understanding, and is not appli- 



cable for the systems confined to the several molecular 
diametersi2iS£ typical for AFM tip-substrate separation. 
Arguably, the approaches based on the Kelvin equation 
do not take into account strongly adsorbed layers, and 
neglect the density oscillations near confining surfaces^ 
There is no molecular-level model describing the effects 
of external electric field to the best of our knowledge. 

The goal of this work is to develop a versatile de- 
scription of the water condensation in proximity of a 
nanoscale asperity and investigate the influence of spa- 
tially non-uniform external electric field, with a focus 
on AFM nanolithography. Here we concern with the 
case of isotropic fluid confined in asymmetric nano junc- 
tion in the absence of the electric field to develop a ba- 
sic framework for the case of field-induced orientational 
anisotropy. 

Typical molecular approaches, applicable to the prob- 
lem of our interest, include computer simulations of the 
Monte Carlo^i and molecular dynamics^ types, and a 
density functional theory (DFT)i^*24 A grand canonical 
Monte Carlo simulation of water condensation in the con- 
text of the dip-pen nanolithography has been reported in 
a series of works by Jang et.ali&lZ^ including clarifi- 
cation of fluctuations influence on the meniscus width. 
The effects of relative humidity, AFM tip curvative, tip 
and substrate wetting properties on the meniscus size and 
capillary force have been studied. 

Unlike other studies, our study concerns with DFT ap- 
proach to model real systems in which the fluctuations 
of meniscus's width, as established in^& are not signifi- 
cant. A mean field character of the DFT is appropriate 
since the fluid is considered far from the critical point. 
Menisci instabilities related to the tip-substrate distance 
variations dependent on the AFM lever thermal fluctu- 
ations, and on other factors could be described in the 
framework of DFT formulation, however, it will be the 
subject of a separate study focused on the specifics of 
AFM-assisted nanolithography. 

Another reason for DFT choice with respect to Monte 
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Carlo technique is the convenience of generalization to 
the case of the electric field induced anisotropy. The 
Monte Carlo simulations with additional orientational 
degrees of freedom would require special cautions to avoid 
local minima and to ensure proper sampling of the com- 
plete configurational space, along with increased com- 
putational time, especially in the cases of the low ac- 
ceptance rates for the regions filled with a dense liquid. 
Molecular dynamics seems to be less convenient as well, 
because of the variable number of particles dictated by 
the grand canonical type of the problem, and low densi- 
ties in the vapor phase potentially providing insufficient 
number of collisions for proper statistical averaging. Ad- 
ditionally, formulation of the analytical expressions, pos- 
sible in the framework of DFT, is highly desirable. Fi- 
nally, another advantage of DFT is a consistency in map- 
ping of the 3D geometry into the quasi-2D description. 
This avoids errors in quantitative prediction of the solva- 
tion forces, present in the case of 2D lattice Monte Carlo 
simulation^ 

The paper is organized as follows. Section [H] is en- 
titled "Model and Methodology". The first part of the 
section describes the model including non-local density 
functional formulation and specifics of fluid, while the 
second part outlines the computational issues related to 
the form of the integral equation for density distribu- 
tion. We discuss multidimensional integrals handling, 
and describe the iterative procedure developed. The Sec- 
tion llll Al concerns with elucidation of the model parame- 
ters of the fluid from the bulk behavior. The Section llll Bl 
presents the results of calculations. The meniscus's size 
and its variation with respect to the relative humidity 
and wall separation are discussed. The Section Till Cl ad- 
dresses determination of the fluid- wall interaction param- 
eters to link the model's calculations to particular mate- 
rial systems. The procedure is based on the polymolec- 
ular adsorption isotherms. The Section ITVl presents the 
summary. 



FIG. 1: Schematic presentation of the nanoscale junction con- 
sisting of a planar surface and a spherical asperity with cur- 
vature radius R, separated by the distance t. The fluid is 
nonuniformly distributed in the junction. The model resem- 
bles AFM tip dwelling above a sample surface. 



The intrinsic Helmholtz free energy functional, F [p (r)], 
incorporates the properties of fluid. In general case, the 
functional depends on the external potential V ex t (r) if 
the intermolecular forces and/or intramolecular degrees 
of freedom are affected by the external field. However, we 
defer this discussion to the second part of the study, and 
consider F [p (r)] invariant with respect to V ex t (r). The 
external potential describes the solid surfaces in contact 
with fluid, and induces inhomogeneiety in the system. 
The chemical potential p is related to the ambient hu- 
midity. 

The Helmholtz free energy functional F [p (r)] can be 
presented as 

F [p {r)\ = F id [p (r)] + AF rep [p (r)] + AF attr [p (r)] , (2) 

where F id [p (r)] corresponds to the ideal system without 
intermolecular interactions, and the other two terms de- 
scribe additional intermolecular repulsion and attraction, 
respectively. The ideal part of the Helmholtz free energy 
functional is given explicitly by 



II. MODEL AND METHODOLOGY 

A. Model description 

The nano-scale junction consists of planar surface and 
spherical asperity of the radius R, separated from the 
surface at the distance t (Fig. The choice of the 
cylindrical coordinate system is related to the azimuthal 
symmetry. The fluid is nonuniformly distributed in the 
junction. 

The purpose of the equilibrium DFT applied to the 
inhomogeneous fluid is to find the spatial distribution for 
the number density of molecules, p (r), by minimizing the 
grand potential f2, presented as a functional of p (r) for 
the grand canonical (T, V, p) ensembl o 23 i 24 : 



F ld [p(r)] = kT f p{r) (in (\ 3 th p(r)) 



f 



^) dr, 
kT ) 



(3) 

Xth is a thermal wavelength, k is the Boltzmann constant, 
T is the temperature and the term pi nt {kT) 1 has been 
added to incorporate internal degrees of freedom of the 
water molecules. 

The repulsive free energy functional is expressed as 



AF rep [p (r)} = kT p{r) Aif> ha (p (r)) dr. 



(4) 



Here Aiphs is an excess free energy density for the repul- 
sive fluid evaluated at every spatial point as a function of 
the smoothed (coarse-grained) density (p(r)). The latter 
is an average of the local density p (r) , taken over a small 
domain using the weighting function lu: 



« [p (r)} = F[p (r)} + / p (r) (V ext (r) - p) dr (1) 



p{f) 



p\r 



dr' 



(5) 
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The weighting function may depend not only on the spa- 
tial coordinate but also be a functional of p, depend- 
ing on the particular non-local DFT formulation^ We 
use the repulsive intermolecular potential with the hard 
sphere diameter a and the Carnahan-Starling free energy 
density 



4i] - 3ij 2 



(6) 



where rj is the packing fraction of the hard spheres cor- 
responding to the number density p: 



V 



(J) 



For the sake of simplicity we choose a generalized van der 
Waals model 23 ' 26 with 



(\ f - ? \) = 4^ 9 (' 7 -\"- p \) 



(8) 



where 9 is a Heaviside step function. While the more so- 
phisticated approaches, such as the Tarazona's model^! 
have better accuracy for describing density profiles of the 
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hard sphere systems, the choice of the generalized van der 
Waals model is justified in the systems with strong inter- 
molecular attraction^ which is the case of our interest. 

In the framework of the mean field approximation, the 
third term in Eq. describing an attractive part of the 
free energy functional can be written as: 



drdr 1 



AF attr [p(r)\ 

"ill \ ' • 

where <& a ttr {f) is the attractive pair intermolecular po- 
tential. We consider the water-like fluid for which 
&attr {r) is a square well potential characterized by two 
effective parameters such as the depth (e) , and the width 
(d) of the well. This choice is related to the fact that 
in water, the short-range hydrogen bonding dominates 
the longer-range dipole-dipole and higher order interac- 
tions. Similar forms of attractive potentials were used 
beforei 28 i 29 

The minimization of the functional given by Eq. 
with F [p (f)] presented above by Eqs. ©-© leads to 
the following integral equation for the equilibrium density 
profile p (r) : 



kT (In (\l h p (r)) - 1) + p, int + &iphsP (f) + As [p(r),p (r)] + J ln t [p (r)] + V ext (r) = p, 



(10a) 



Jhs [p (r) , P (r)} 



Ana 3 



dr' 



(10b) 



Jint [p(?)]=-ej 9 (\f~ r'\ -v)e(d- r-r' ),,[,■') <lr'. 

I 



(10c) 



where a prime at Aiphs in <|10bfl denotes the differentia- 
tion with respect to p. The Eq. (|10a|l is to be solved for 
0<r<L,0<z<t + R- y/R 2 - r 2 , and < <j) < 2ir. 
The boundary r — Lis chosen such that the region within 
which the liquid bridge may be formed corresponds to 
r < L. 

An explicit form of the external potential, V ext (f) de- 
scribing interaction of the fluid with the confining sur- 
faces is the final piece of the information required for 
closure of the model's description. For that we use a 
two-parametric exponential attraction potential on top 
of the hard wall repulsion: 



V ext (x) = 



-e s exp(- 

DO, 



ax 



x>0 
x < 



(11) 



where x is the distance from the surface. The parameter 
e s is the work required to displace a water molecule from 
x = to infinity corresponding to the binding energy of 



water at given surface. 



B. Computational part 

The triple integrals l|10b|) and l|10cfl must be evaluated 
over 3-D volume. Fortunately, the character of their an- 
gular dependence in a step- function form allows reducing 
them to the double integrals below: 



u (r) 6 (p — \f— ro|) dr 



rdr J dz ■ 2Xu (r, z) , 
(12a) 



arccos (A) , -1 < A < 1 
0, A>1 
7T, A < -1 



(12b) 



4 



III. RESULTS AND DISCUSSION 



A 



(z - z Q ) 2 + r 2 + rl 
2rr n 



(12c) 



The iterative method used to solve nonlinear integral 
Eq. HI U| l deserves special description. Simple successive 
approximation procedure normally fails for the equations 
similar to I|1U|) even for one-dimensional problems as the 
computation diverges after the few iterationsi^SiSi Sev- 
eral approaches were used to overcome this difficulty for 
the one-dimensional case. Tarazona££ used the mixing 
of the new and previous iterative values with a space- 
dependent mixing function. While this approach fixes 
divergence, it requires a careful selection of the mixing 
function through trial end error, and converges slowly. 
Vanderlick et.alm^ used a uniform discretization of the 
domain and converted an integral equation into the sys- 
tem of nonlinear algebraic equations. This works well for 
one-dimensional problems, but is subjective to the choice 
of the discretization method. 

In our model, the spatial domain for which Eq. (|1(J|) has 
to be solved is two dimensional and asymmetric, which 
requires a general approach. Below we suggest an iter- 
ative scheme conceptually similar to Newton's iterative 
method for the systems of algebraic equations. The gen- 
eral form of Eq. (|10fl can be written as: 



lnp + $(/?) = B, 



(13) 



where $ (/?) has a complicated dependence with respect 
to p. The iterative scheme is set up in the form 



In 



d<f> 

~d~ P 



B. 



(14) 



Here is an approximation for p (r, z) obtained at k th 
iteration. The integral terms appearing in $ and its 
derivative are calculated simultaneously. The equation 
to be solved at m th iteration has the following form: 



hV m+1 > +Dp (m+1 ^ = C, 



(15) 



A. Bulk fluid and model parameters 

For the case of the uniform bulk fluid, Eq. (|10fl reduces 
to the following form: 

/i = kT (In (\ 3 th p) - 1) + fi int + A^ hs + pA<P' hs 




To determine the parameters <r, d and e in Eq. (|17fl we 
use the characteristics of water at the critical point de- 
termined by the following thermodynamic conditions 30 : 



Oft 
dp 



0, 



d 2 p 
dp 2 



(18) 



The derivatives are evaluated at the critical temperature 
T c (647 K for water). The solution of these equations 
with /i explicitly given by Eq. I|17fl determines two groups 
of the parameters: 



cr 3 p c = 0.2372, 





10.3882. (19) 



The density of water under critical point conditions, p c , 
results in a — 2.8A, which is commonly taken as a hard 
sphere diameter of water molecule^ Separation of the 
parameters e and d in the second group is not necessary 
as they are always grouped together in this form. 

The second consideration for the bulk fluid deals with 
liquid-vapor coexistence at the temperature T = 300 K. 
Two coexistence equations are derived assuming the 
chemical potential and pressure equal for both phases. 
The chemical potential is given by Eq. I|17fl . and the pres- 
sure is presented as a function of number density, P (p) , 
using the expression for the grand canonical potential, Q, 
for the bulk fluid: 



(7° 

kf 



/>-^!Q + ln(a 3 p)+^ 




1 



kT 



(20) 



and its solution can be presented via the Lambert W- 
function: 



exp(C- W(De c )) 



(16) 



We found that the iterative procedure given above 
through Eqs. i|14[l - i|16|) i s ver y efficient so that the conver- 
gence is often accomplished in less than 10-20 iterations. 



where dimensionless parameter Q was introduced as 

Q = 31n^)-1 + ^. (21) 

After solving the equations above, we find <r 3 p v = 1.86 x 
10~ 4 and Q — 5.6 selecting the density of water in the 
liquid phase at 300 K, a 3 pi — 0.731. The vapor phase 
density p v is used as a reference value to calculate chem- 
ical potential p. for a given relative humidity H r . 
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FIG. 2: Water condensation between planar surface and 
spherical asperity of the radius R = 36a (10 nm). Shaded 
area corresponds to the formation of liquid meniscus. No liq- 
uid bridge is formed for given H r when t reaches ordinate 
of the open circle. The parameters of fluid interaction with 
confining walls are e s e _1 = 1 and aa = 2. 



B. Liquid meniscus condensation in the nanoscale 
junction 

The fluid density distribution was calculated, using the 
methodology above, in the geometry illustrated in Fig.^ 
for different values of relative humidity H r and separa- 
tion between confining walls t, at T = 300 K. It was 
found that for certain values of H r and t, the fluid tends 
to condense in the region near r = 0, forming a liquid 
meniscus that connects the confining surfaces. The ten- 
dency for liquid meniscus condensation increases with H r 
and diminishes with the growth of t. The region in H r — t 
coordinates, corresponding to the meniscus formation, is 
depicted in Fig. [21 

Fig. 13 presents typical density profiles: a 3 p as a func- 
tion of the axial coordinate za~ x , for different radial po- 
sitions r. The density profile for a given r is mainly 
determined by the separation between the walls at that 
point, providing that the fluid- wall interaction parame- 
ters are fixed. Inside the meniscus, the density is higher 
near the walls and exhibits maxima in-between (Fig.[3J). 
Significant density redistribution occurs near the menis- 
cus boundary. For the longer radial distances, the liquid 
wets confining surfaces according to the specified values 
of e s and a. 

An essential question is how to define the size of liquid 
meniscus based on the computed fluid density distribu- 
tion, or what factors define liquid/ vapor boundary for a 
smooth density distribution in the asymmetric geometry. 
The answer can be found should we specify the phenom- 
ena of interest and its description based on the dimen- 
sions of the liquid bridge. The description of the charge 
transport in nanoscale junction and calculation of the 
forces exerted by fluid on the confining walls may require 






FIG. 3: Fluid density profiles between planar surface and 
spherical asperity of the radius R — 36a (10 nm), separated by 
t = 5.5<r, for three different values of the radial coordinate r. 
Relative humidity is 90%, and parameters of fluid interaction 
with solid walls are e s e _1 = 1 and aa — 2. 



the usage of different effective meniscus radii. A certain 
approximation is required since this work studies general 
peculiarities of meniscus formation under variable condi- 
tions. In earlier work related to meniscus formation^ an 
ad hoc chosen threshold of half filled lattice site has been 
used. In this work, we use a single threshold value of 
the number density to separate liquid and vapor regions, 
and define the boundary. We assume that if the interac- 
tion energy of water molecule located at a given point of 
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FIG. 4: Water meniscus formed between planar surface and 
spherical asperity of the radius R = 36a (10 nm), separated 
by t = 5.5cr. Relative humidity is 90%, and parameters of 
fluid interaction with solid walls are e s e _1 = 1 and aa — 2. 



space with its neighbors is equal or less than kT, then it 
belongs to vapor phase. The corresponding quantitative 
criterion can be presented using Eq. (|10c|) as 

-Jim (p) = kT, (22) 

or 

This results in the threshold value between the liquid 
and vapor phases to be a 3 pth = 0.16. For comparison, 
we found that the bulk liquid density was a 3 pi — 0.731, 
while the vapor density was a 3 p v — 1.86 x 10~ 4 (as de- 
scribed in Section UlI Al above) . 

Fig. 21 presents the liquid meniscus that corresponds to 
the density profiles plotted in Fig. |3| The spaces on the 
plot filled with dots correspond to the regions (domain 
discretization points) for which the fluid density exceeds 
the threshold value p t h defined above. 

Fig. illustrates the variation of the meniscus's size 
with respect to the relative humidity H r . Again, the 
spaces on the plots filled with dots correspond to the re- 
gions for which the fluid density exceeds pth- The menis- 
cus boundary shifts to larger values of r as humidity 
grows. 

Fig. shows similar plots for the variable separation t 
between the walls. The meniscus narrows (as can be seen 
for ta^ 1 =2 — 3) and eventually disappears (ta^ 1 = 3.5) 
as t increases. 
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FIG. 5: Variation of the meniscus size with relative humidity 
H r . Filled space on the plots corresponds to the regions where 
fluid density exceeds the threshold value of a 3 pth = 0.16. 
Other relevant parameters are: t = 3<r, e s e _1 = 1 and aa — 2. 
Curvature radius of the asperity is 7? = 36a. 



C. Fluid - wall interaction parameters: An 
Estimation 

The parameters, e s and a (Eq. I|ll|l ) characterizing the 
energy of the fluid interaction with solid surfaces (e s ), 
and the spatial range (a) of the corresponding poten- 
tial, must be estimated for specific surfaces to apply the 
results of the calculations to real physical systems. We 
suggest using the polymolecular adsorption isotherms for 
this purpose since the data of water adsorption is more 
abundant than those related to other measurements. 
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FIG. 7: Calculated adsorption isotherms for the model fluid 
on surfaces with different values of e a e~ : 1 (dash-dot line), 
1.5 (dotted line), 2 (solid line), and with aa — 2. V is num- 
ber of molecules adsorbed per unit surface area. As H r ap- 
proaches 100%, the onset of (integrable) singularity is ob- 
served (not shown). 
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FIG. 6: Variation of the meniscus size with separation t be- 
tween confining surfaces. Filled space on the plots corre- 
sponds to the regions where fluid density exceeds the thresh- 
old value of a 3 pth = 0.16. Relative humidity is 50%, and 
parameters of fluid interaction with solid walls are e s e _1 = 1 
and aa — 2. Curvature radius of the asperity is R — 36a. 



The fluid-wall interaction potential is two- 

parametric. First, we restrict the potential to act in the 
region of thickness about a near the surface by setting 
aa = 2. Then the adsorption isotherms for the model 
fluid on surfaces with different values of e s are calculated. 

The number of molecules adsorbed per surface area of 
a 2 (integrated or Gibbs adsorption) with respect to the 
relative humidity H r for T = 300K is presented in Fig. [7| 
A singularity is observed (not shown in the figure) as the 
humidity (H r ) approaches 100%, consistent with thermo- 
dynamical requirement The isotherms exhibit typical 
step-like structure {est" 1 ~ 1) for the intermediate hu- 
midity values, similar to the results of calculations based 
on three-dimensional Ono-Kondo lattice modei^ and ex- 
perimental data for multilayer adsorption of water ^ The 
parameter e s for specific surface of interest may thus be 
evaluated by fitting the adsorption isotherm. 

The above considerations imply both geometric and 
energetic uniformity of the surface. The geomet- 
ric one means the average surface roughness is much 
smaller than the separation distance (t) between sur- 
faces. The latter assumption can be adapted to ener- 
getically nonuniform surfaces by taking the effective (av- 
erage) value of e s and requiring that the radial region of 
interest is large enough to sample most of the adsorp- 
tion energy distribution. Additionally, it is assumed that 
no significant dissociation of water occurs on the surface, 
and adsorption is reversible. 
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IV. SUMMARY 

A spatial distribution of the water-like fluid in the 
asymmetric nanoscale junction formed in proximity of 
a nanoscale asperity has been studied using an approach 
based on non-local density functional theory. The ap- 
proach presents a closed set of procedures modeling the 
fluid resembling water. The computational part of the 
approach is adapted for asymmetric geometry of the 
model comprising two solid surfaces. The basic frame- 
work for the density functional description of the case 
with electric field induced anisotropy in the fluid has been 
developed. The route for quantitative predictions on the 
water meniscus formation near atomic force microscope 
tip dwelling above a sample surface is one of the impor- 
tant applications of the developed methodology. 

Our results suggest further development of the model- 
ing in three directions. The first direction would be a self- 



consistent description for electrically charged or biased 
asperity comprising electric field and fluid density distri- 
butions. The second direction would be the calculation 
of forces acting on the nano-asperity based on the density 
functional formalism, which would have applications in 
the context of AFM-assisted nanolithography. The third 
direction would concentrate on the study of dynamical 
phenomena at the nano-asperity including the kinetics of 
capillary condensation under highly nonuniform external 
potentials. 
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